Acanthophora Pmax & NPQ Model Script

Acanthophora spicifera PHOTOSYNTHESIS (Pmax and NPQ) Analysis, Script Chunks, and Plots

This is the analysis of two runs of Acanthophora spicifera salinity and nutrient experiments conducted on the lanai in St. John 616 in November and December 2023. These experiments incorporated two salinity/nutrient levels and two temperature levels. Pmax, maximum NPQ and delta NPQ are used for these analyses. The max NPQ used is that which also satisfies the need for yield to be greater than 0.1. deltaNPQ takes that maximum value and subtracts the initial value of NPQ that is recorded at the second sat pulse when Epar is at 66.

Load the appropriate packages

library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
#library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load this file for normalized to quantum efficiency of photosynthesis per Silsbe and Kromkamp

acan_ps <- read.csv("../transformed/acan_ek_alpha_normalized_2024.csv")

Assign run as a factor

acan_ps$run <- as.factor(acan_ps$run)

Assign temperature as a factor

acan_ps$temp <- as.factor(acan_ps$temp)

Assigns treatment as characters from integers then to factors

acan_ps$treatment <- as.factor(as.character(acan_ps$treatment))

Assign deltaNPQ as numerical

acan_ps$deltaNPQ <- as.numeric(acan_ps$deltaNPQ)

Use Day 9 for final analysis

acan_day9 <- subset(acan_ps, rlc_day == 9)
acan_day9$treatment_graph[acan_day9$treatment == 0] <- "1) 35ppt/0.5umol"
acan_day9$treatment_graph[acan_day9$treatment == 1] <- "2) 28ppt/14umol"

acan_day9$temperature_graph[acan_day9$temp == 24] <- "24°C"
acan_day9$temperature_graph[acan_day9$temp == 28] <- "28°C"

##Pmax Plot histogram for Pmax

acan_day9 %>% ggplot(aes(pmax)) +
  geom_histogram(binwidth=5, fill = "goldenrod1", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

Run model

acan_day9_pmax_model <- lmer(formula = pmax ~ treatment + temp + (1 | run) + (1 | plantID) + (1 | rlc_order), data = acan_day9)

Construct null model to perform likelihood ratio test REML must be FALSE

acan_pmax_treatment_null <- lmer(formula = pmax ~ temp + (1 | run) + (1 | plantID) + (1 | rlc_order), data = acan_day9, REML = FALSE)
acan_pmax_model2 <- lmer(formula = pmax ~ treatment + temp + (1 | run) + (1 | plantID) + (1 | rlc_order), data = acan_day9, REML = FALSE)
anova(acan_pmax_treatment_null, acan_pmax_model2)

## Data: acan_day9
## Models:
## acan_pmax_treatment_null: pmax ~ temp + (1 | run) + (1 | plantID) + (1 | rlc_order)
## acan_pmax_model2: pmax ~ treatment + temp + (1 | run) + (1 | plantID) + (1 | rlc_order)
##                          npar    AIC    BIC  logLik deviance  Chisq Df
## acan_pmax_treatment_null    6 729.65 745.03 -358.82   717.65          
## acan_pmax_model2            7 731.48 749.43 -358.74   717.48 0.1696  1
##                          Pr(>Chisq)
## acan_pmax_treatment_null           
## acan_pmax_model2             0.6805

acan_pmax_temperature_null <- lmer(formula = pmax ~ treatment + (1 | run) + (1 | plantID) + (1 | rlc_order), data = acan_day9, REML = FALSE)
acan_pmax_model3 <- lmer(formula = pmax ~ treatment + temp + (1 | run) + (1 | plantID) + (1 | rlc_order), data = acan_day9, REML = FALSE)
anova(acan_pmax_temperature_null, acan_pmax_model3)

## Data: acan_day9
## Models:
## acan_pmax_temperature_null: pmax ~ treatment + (1 | run) + (1 | plantID) + (1 | rlc_order)
## acan_pmax_model3: pmax ~ treatment + temp + (1 | run) + (1 | plantID) + (1 | rlc_order)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## acan_pmax_temperature_null    6 742.95 758.33 -365.47   730.95          
## acan_pmax_model3              7 731.48 749.43 -358.74   717.48 13.467  1
##                            Pr(>Chisq)    
## acan_pmax_temperature_null               
## acan_pmax_model3            0.0002428 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Make residual plots of the data for Acanthophora

hist(resid(acan_day9_pmax_model))

plot(resid(acan_day9_pmax_model) ~ fitted(acan_day9_pmax_model))

qqnorm(resid(acan_day9_pmax_model))
qqline(resid(acan_day9_pmax_model))

Check the performance of the model

performance::check_model(acan_day9_pmax_model)

Get R squared

r.squaredGLMM(acan_day9_pmax_model)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##            R2m       R2c
## [1,] 0.2337334 0.5828728

Make plots and tables for the data

tab_model(acan_day9_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	pmax
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	51.04	2.95	45.19 – 56.90	17.32	<0.001	89.00
treatment [1]	-1.04	2.71	-6.43 – 4.35	-0.38	0.703	89.00
temp [28]	-11.99	3.04	-18.03 – -5.96	-3.95	<0.001	89.00
Random Effects
σ²	65.35
τ₀₀ _plantID	22.46
τ₀₀ _{rlc_order}	27.78
τ₀₀ _run	4.46
ICC	0.46
N _run	2
N _plantID	48
N _{rlc_order}	24
Observations	96
Marginal R² / Conditional R²	0.234 / 0.583

plot(allEffects(acan_day9_pmax_model))

acan_day9 %>% ggplot(aes(treatment_graph, pmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = temp), position = "jitter", show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "28ppt/14umolN")) + 
  ylim(-1, 90) + stat_mean() + 
  scale_color_manual(values = c("brown", "brown2")) +
  geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

acan_day9 %>% ggplot(aes(temperature_graph, pmax)) +
  geom_boxplot(size=0.5) +
  geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) +
  labs(x="temperature", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Acanthophora spicifera") +
  ylim(-1, 90) + stat_mean() + 
  scale_color_manual(values = c("brown4", "brown1")) +
  geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means for pmax

acan_day9 %>% group_by(treatment) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 2 × 2
##   treatment  mean
##   <fct>     <dbl>
## 1 0          45.0
## 2 1          44.0

acan_day9 %>% group_by(temp) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 2 × 2
##   temp   mean
##   <fct> <dbl>
## 1 24     50.5
## 2 28     38.5

##Linear regression Pmax vs Growth——————————————————————————–

Add growth rate from other dataset to this one and subset by species

growth_rate <- read.csv("../input/acanthophora_growth_r1.csv")
growth_rate$treatment <- as.factor(growth_rate$treatment)

Make a new column for weight change (difference final from initial)

growth_rate$growth_rate_percent <- (growth_rate$final_weight - growth_rate$initial_weight) / growth_rate$initial_weight * 100
acan_day9$growth_rate_percent <- growth_rate$growth_rate_percent
acan_day9$secondary_apices <- growth_rate$secondary_apices
acan_day9_sub <- subset(acan_day9, secondary_apices != "200")

Plot a regression between the photosynthetic independent variables of interest and growth rate

acan_growth_pmax_graph <- ggplot(acan_day9_sub, aes(x=pmax, y=growth_rate_percent)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Acanthophora spicifera Pmax vs 9-Day Growth (%)", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "9-Day Growth (%)") + stat_regline_equation(label.x = 10, label.y = 65) + stat_cor(label.x = 10, label.y = 60)
acan_growth_pmax_graph

## `geom_smooth()` using formula = 'y ~ x'

Plot a regression between the photosynthetic independent variables of interest and apices/axis

acan_apices_pmax_graph <- ggplot(acan_day9_sub, aes(x=pmax, y=secondary_apices)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Acanthophora spicifera Pmax vs # Secondary Apices", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "Secondary Apices (#)") + stat_regline_equation(label.x = 10, label.y = 65) + stat_cor(label.x = 10, label.y = 60)
acan_apices_pmax_graph

## `geom_smooth()` using formula = 'y ~ x'

NPQ analysis using max values that also satisfy the yield>0.1 argument

Plot data histogram

acan_day9 %>% ggplot(aes(maxNPQ_Ypoint1)) +
  geom_histogram(binwidth=0.05, fill = "violet", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

Run model without interaction between the treatments and temperature

acan_day9_maxNPQ_model <- lmer(formula = maxNPQ_Ypoint1 ~ treatment + temp + (1 | run) + (1 | plantID), data = acan_day9)

## boundary (singular) fit: see help('isSingular')

Construct null model to perform likelihood ratio test REML must be FALSE

acan_maxNPQ_treatment_null <- lmer(formula = maxNPQ_Ypoint1 ~ temp + (1 | run) + (1 | plantID), data = acan_day9, REML = FALSE)

## boundary (singular) fit: see help('isSingular')

acan_maxNPQ_model2 <- lmer(formula = maxNPQ_Ypoint1 ~ treatment + temp + (1 | run) + (1 | plantID), data = acan_day9, REML = FALSE)

## boundary (singular) fit: see help('isSingular')

anova(acan_maxNPQ_treatment_null, acan_maxNPQ_model2)

## Data: acan_day9
## Models:
## acan_maxNPQ_treatment_null: maxNPQ_Ypoint1 ~ temp + (1 | run) + (1 | plantID)
## acan_maxNPQ_model2: maxNPQ_Ypoint1 ~ treatment + temp + (1 | run) + (1 | plantID)
##                            npar     AIC     BIC logLik deviance  Chisq Df
## acan_maxNPQ_treatment_null    5 -73.970 -61.148 41.985   -83.97          
## acan_maxNPQ_model2            6 -88.553 -73.167 50.277  -100.55 16.584  1
##                            Pr(>Chisq)    
## acan_maxNPQ_treatment_null               
## acan_maxNPQ_model2          4.655e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

acan_maxNPQ_temperature_null <- lmer(formula = maxNPQ_Ypoint1 ~ treatment + (1 | run) + (1 | plantID), data = acan_day9, REML = FALSE)

## boundary (singular) fit: see help('isSingular')

acan_maxNPQ_model3 <- lmer(formula = maxNPQ_Ypoint1 ~ treatment + temp + (1 | run) + (1 | plantID), data = acan_day9, REML = FALSE)

## boundary (singular) fit: see help('isSingular')

anova(acan_maxNPQ_temperature_null, acan_maxNPQ_model3)

## Data: acan_day9
## Models:
## acan_maxNPQ_temperature_null: maxNPQ_Ypoint1 ~ treatment + (1 | run) + (1 | plantID)
## acan_maxNPQ_model3: maxNPQ_Ypoint1 ~ treatment + temp + (1 | run) + (1 | plantID)
##                              npar     AIC     BIC logLik deviance  Chisq Df
## acan_maxNPQ_temperature_null    5 -90.469 -77.647 50.234  -100.47          
## acan_maxNPQ_model3              6 -88.553 -73.167 50.277  -100.55 0.0847  1
##                              Pr(>Chisq)
## acan_maxNPQ_temperature_null           
## acan_maxNPQ_model3               0.7711

Make residual plots of the data for Acanthophora

hist(resid(acan_day9_maxNPQ_model))

plot(resid(acan_day9_maxNPQ_model) ~ fitted(acan_day9_maxNPQ_model))

qqnorm(resid(acan_day9_maxNPQ_model))
qqline(resid(acan_day9_maxNPQ_model))

Check the performance of the model

performance ::check_model(acan_day9_maxNPQ_model)

Get R-squared

r.squaredGLMM(acan_day9_maxNPQ_model)

##            R2m       R2c
## [1,] 0.1508684 0.2570364

Make plots and tables for the data

tab_model(acan_day9_maxNPQ_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	maxNPQ_Ypoint1
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	0.51	0.03	0.46 – 0.56	19.30	<0.001	90.00
treatment [1]	-0.12	0.03	-0.18 – -0.07	-4.38	<0.001	90.00
temp [28]	-0.01	0.03	-0.07 – 0.05	-0.28	0.776	90.00
Random Effects
σ²	0.02
τ₀₀ _plantID	0.00
τ₀₀ _run	0.00
ICC	0.13
N _run	2
N _plantID	48
Observations	96
Marginal R² / Conditional R²	0.151 / 0.257

plot(allEffects(acan_day9_maxNPQ_model))

acan_day9 %>% ggplot(aes(treatment_graph, maxNPQ_Ypoint1)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = temp), position = "jitter", show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Day 9 NPQmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "28ppt/14umolN")) + 
  ylim(0, 1) + stat_mean() + 
  scale_color_manual(values = c("brown", "brown2")) +
  geom_hline(yintercept=0.5, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

acan_day9 %>% ggplot(aes(temperature_graph, maxNPQ_Ypoint1)) +
  geom_boxplot(size=0.5) +
  geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) +
  labs(x="temperature", y= "Day 9 NPQmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Acanthophora spicifera") +
  ylim(0, 1) + stat_mean() + 
  scale_color_manual(values = c("olivedrab3", "maroon4")) +
  geom_hline(yintercept=0.5, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means for NPQ

acan_day9 %>% group_by(treatment) %>% summarise_at(vars(maxNPQ_Ypoint1), list(mean = mean))

## # A tibble: 2 × 2
##   treatment  mean
##   <fct>     <dbl>
## 1 0         0.505
## 2 1         0.382

acan_day9 %>% group_by(temp) %>% summarise_at(vars(maxNPQ_Ypoint1), list(mean = mean))

## # A tibble: 2 × 2
##   temp   mean
##   <fct> <dbl>
## 1 24    0.448
## 2 28    0.439

lot a regression between the photosynthetic independent variables of interest and growth rate

acan_growth_maxNPQ_graph <- ggplot(acan_day9_sub, aes(x=maxNPQ_Ypoint1, y=growth_rate_percent)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Acanthophora spicifera NPQmax vs 9-Day Growth (%)", x = "NPQmax (rel. units)", 
       y = "9-Day Growth (%)") + stat_regline_equation(label.x = 0.25, label.y = 38) + stat_cor(label.x = 0.25, label.y = 40)
acan_growth_maxNPQ_graph

## `geom_smooth()` using formula = 'y ~ x'

Plot a regression between the photosynthetic independent variables of interest and apices/axis

acan_apices_NPQmax_graph <- ggplot(acan_day9_sub, aes(x=maxNPQ_Ypoint1, y=secondary_apices)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Acanthophora spicifera NPQmax vs # Secondary Apices", x = "NPQmax (rel. units)", 
       y = "Secondary Apices (#)") + stat_regline_equation(label.x = 0.25, label.y = 125) + stat_cor(label.x = 0.25, label.y = 128)
acan_apices_NPQmax_graph

## `geom_smooth()` using formula = 'y ~ x'

Delta NPQ

Plot histogram

acan_day9 %>% ggplot(aes(deltaNPQ)) +
  geom_histogram(binwidth=.05, fill = "violet", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

Run model without interaction between the treatments and temperature

acan_day9_deltaNPQ_model <- lmer(formula = deltaNPQ ~ treatment + temp + (1 | plantID), data = acan_day9)

Construct null model to perform likelihood ratio test REML must be FALSE

acan_deltaNPQ_treatment_null <- lmer(formula = deltaNPQ ~ temp + (1 | plantID), data = acan_day9, REML = FALSE)

## boundary (singular) fit: see help('isSingular')

acan_deltaNPQ_model2 <- lmer(formula = deltaNPQ ~ treatment + temp + (1 | plantID), data = acan_day9, REML = FALSE)
anova(acan_deltaNPQ_treatment_null, acan_deltaNPQ_model2)

## Data: acan_day9
## Models:
## acan_deltaNPQ_treatment_null: deltaNPQ ~ temp + (1 | plantID)
## acan_deltaNPQ_model2: deltaNPQ ~ treatment + temp + (1 | plantID)
##                              npar     AIC     BIC logLik deviance  Chisq Df
## acan_deltaNPQ_treatment_null    4 -163.33 -153.07 85.665  -171.33          
## acan_deltaNPQ_model2            5 -180.87 -168.04 95.433  -190.87 19.537  1
##                              Pr(>Chisq)    
## acan_deltaNPQ_treatment_null               
## acan_deltaNPQ_model2          9.868e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

acan_deltaNPQ_temperature_null <- lmer(formula = deltaNPQ ~ treatment + (1 | plantID), data = acan_day9, REML = FALSE)
acan_deltaNPQ_model3 <- lmer(formula = deltaNPQ ~ treatment + temp + (1 | plantID), data = acan_day9, REML = FALSE)
anova(acan_deltaNPQ_temperature_null, acan_deltaNPQ_model3)

## Data: acan_day9
## Models:
## acan_deltaNPQ_temperature_null: deltaNPQ ~ treatment + (1 | plantID)
## acan_deltaNPQ_model3: deltaNPQ ~ treatment + temp + (1 | plantID)
##                                npar     AIC     BIC logLik deviance Chisq Df
## acan_deltaNPQ_temperature_null    4 -180.85 -170.59 94.426  -188.85         
## acan_deltaNPQ_model3              5 -180.87 -168.04 95.433  -190.87 2.013  1
##                                Pr(>Chisq)
## acan_deltaNPQ_temperature_null           
## acan_deltaNPQ_model3                0.156

Make residual plots of the data for Acanthophora

hist(resid(acan_day9_deltaNPQ_model))

plot(resid(acan_day9_deltaNPQ_model) ~ fitted(acan_day9_deltaNPQ_model))

qqnorm(resid(acan_day9_deltaNPQ_model))
qqline(resid(acan_day9_deltaNPQ_model))

Check the performance of the model

performance ::check_model(acan_day9_deltaNPQ_model)

Get R-squared

r.squaredGLMM(acan_day9_deltaNPQ_model)

##            R2m       R2c
## [1,] 0.1847836 0.3356705

Make plots and tables for the data

tab_model(acan_day9_deltaNPQ_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	deltaNPQ
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	0.43	0.02	0.40 – 0.47	26.00	<0.001	91.00
treatment [1]	-0.08	0.02	-0.12 – -0.05	-4.85	<0.001	91.00
temp [28]	-0.03	0.02	-0.07 – 0.01	-1.40	0.164	91.00
Random Effects
σ²	0.01
τ₀₀ _plantID	0.00
ICC	0.19
N _plantID	48
Observations	96
Marginal R² / Conditional R²	0.185 / 0.336

plot(allEffects(acan_day9_deltaNPQ_model))

acan_day9 %>% ggplot(aes(treatment_graph, deltaNPQ)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = temp), position = "jitter", show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Day 9 Delta NPQ (μmols electrons m-2 s-1)", title= "A", subtitle = "Acanthophora spicifera") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "28ppt/14umolN")) + 
  ylim(0, 0.75) + stat_mean() + 
  scale_color_manual(values = c("olivedrab2", "maroon2")) +
  geom_hline(yintercept=0.25, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

acan_day9 %>% ggplot(aes(temperature_graph, deltaNPQ)) +
  geom_boxplot(size=0.5) +
  geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) +
  labs(x="temperature", y= "Day 9 Delta NPQ (μmols electrons m-2 s-1)", title= "A", subtitle = "Acanthophora spicifera") +
  ylim(0, 0.75) + stat_mean() + 
  scale_color_manual(values = c("olivedrab3", "maroon4")) +
  geom_hline(yintercept=0.25, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means for deltaNPQ

acan_day9 %>% group_by(treatment) %>% summarise_at(vars(deltaNPQ), list(mean = mean))

## # A tibble: 2 × 2
##   treatment  mean
##   <fct>     <dbl>
## 1 0         0.420
## 2 1         0.338

acan_day9 %>% group_by(temp) %>% summarise_at(vars(deltaNPQ), list(mean = mean))

## # A tibble: 2 × 2
##   temp   mean
##   <fct> <dbl>
## 1 24    0.394
## 2 28    0.365

Linear regression deltaNPQ vs Growth and Apices

acan_growth_deltaNPQ_graph <- ggplot(acan_day9_sub, aes(x=deltaNPQ, y=growth_rate_percent)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Acanthophora spicifera Delta NPQ vs 9-Day Growth (%)", x = " Delta NPQ (rel. units)", 
       y = "9-Day Growth (%)") + stat_regline_equation(label.x = 0.25, label.y = 38) + stat_cor(label.x = 0.25, label.y = 40)
acan_growth_deltaNPQ_graph

## `geom_smooth()` using formula = 'y ~ x'

Plot a regression between the photosynthetic independent variables of interest and apices/axis

acan_apices_deltaNPQ_graph <- ggplot(acan_day9_sub, aes(x=deltaNPQ, y=secondary_apices)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Acanthophora spicifera Delta NPQ vs # Secondary Apices", x = "Delta NPQ (rel. units)", 
       y = "Secondary Apices (#)") + stat_regline_equation(label.x = 0.25, label.y = 125) + stat_cor(label.x = 0.25, label.y = 128)
acan_apices_deltaNPQ_graph

## `geom_smooth()` using formula = 'y ~ x'

Acanthophora Pmax & NPQ Model Script

Angela Richards Donà

2024-05-15

Acanthophora spicifera PHOTOSYNTHESIS (Pmax and NPQ) Analysis, Script Chunks, and Plots

NPQ analysis using max values that also satisfy the yield>0.1 argument

Delta NPQ